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^ ! ABSTRACT 

^ ' Accretion in protoplanetary discs is thought to be driven by magnetohydrodynamic (MHD) 

^ I turbulence via the magnetorotational instability (MRI). Recent work has shown that a plan- 

■ etesimal swarm embedded in a fully turbulent disc is subject to strong excitation of the veloc- 
' ity dispersion, leading to collisional destruction of bodies with radii Rp < 100 km. Significant 

diffusion of planetesimal semimajor axes also arises, leading to large-scale spreading of the 

■ planetesimal population throughout the inner regions of the protoplanetary disc, in apparent 
p ^ I contradiction of constraints provided by the distribution of asteroids within the asteroid belt, 
[jj . In this paper, we examine the dynamics of planetesimals embedded in vertically stratified 

• turbulent discs, with and without dead zones. Our main aims are to examine the turbulent 

' excitation of the velocity dispersion, and the radial diff'usion, of planetesimals in these discs. 

We employ three dimensional MHD simulations using the shearing box approximation, along 

Q ' with an equilibrium chemistry model that is used to calculate the ionisation fraction of the disc 

( , gas as a function of time and position. Ionisation is assumed to arise because of stellar X-rays, 

C/3 ■ galactic cosmic rays and radioactive nuclei. In agreement with our previous study, we find that 

d ' planetesimals in fully turbulent discs develop large random velocities that will lead to colli- 
sional destruction/erosion for bodies with sizes below 100 km, and undergo radial diffusion 

■ on a scale ~ 2.5 au over a 5 Myr disc life time. But planetesimals in a dead zone experience 
^ I a much reduced excitation of their random velocities, and equilibrium velocity dispersions lie 

. between the disruption thresholds for weak and strong aggregates for sizes Rp < 100 km. We 

OO ' also find that radial diffusion occurs over a much reduced length scale ~ 0.25 au over the disc 

0^ I life time, this being consistent with solar system constraints. We conclude that planetesimal 

CO ■ growth via mutual collisions between smaller bodies cannot occur in a fully turbulent disc. 

' By contrast, a dead zone may provide a safe haven in which km-sized planetesimals can avoid 

, mutual destruction through collisions. 

■ Key words: accretion disks - magnetohydrodynamics (MHD) - methods; numerical - plan- 
• • ' etary systems: formation - planetary systems: protoplanetary disks 



1 INTRODUCTION location in the disc, implying formati on times of =s few x l(f yr at 

„, „ . „ , .... . . „ . „ 5 au. More recent models presented bv lBrauer et al. that ex- 

The formation of planetesimals is a key stage in the formation of . , • , j- • , , • , , ,■ 

amine planetesimal formation in low density turbulent discs, con- 
planetary systems, but at present there is little consensus about , , , ■ , , ■ , j, ■ ■ , , , 

elude that incremental planetesimal lormation is prevented by the 
the processes involved. One class of models envisages an mere- ... . c ^ jiu i »i r 

. . . . " rapid inward migration of metre-sized boulders when they form, 

mental process in which small dust grains collide, stick and set- , ■ , ■ , r • ■ i i i i • i i ■ n- ■ 

combined with fragmentation induced by high velocity collisions 
tie vertically within a protoplanetary disc, with continued growth ,, . . , ■ ,, ,■ 

. . caused by differential radial migration. Large km-sized bodies were 

through coagulation eventually forming km-sized planetesimals i i . ■ j ,.i j , ■_ i- n i i r 

I . " —fr- — — — 1 1 "" ' L i I . " , ii- I PTAnri ivj J 1 '-"^'y obtained in these models when unrealistically large values for 



dWeidenschilling & Cuzzil Il993l : IWeidenschillind llOOOl) . Models u -,, rr jj^ - r 

1 . the critical velocity tor fragmentation were adopted, and a signif- 

discussed by Weidenschillina ( 20001 ) ■ that assume the disc is lami- ■ , , ■ , , , , 



icant increase in the dust-to-gas ratio above solar abundance was 
nar and has a surface density a few times larger than the minimum j a i i i ■ ■_ j i r .i i j . 

u — 1 I 1° assumed. A local particle-in-a-box model oi growth through dust 

mass solar nebula model of IHavashil ( il981f) . predict planetesimal , . , j u Ft TTl STfvuh tu t ■ 



formation times that are a few thousand local orbital periods at each 
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coagulation presented by lZsometaLlbOld) . that incorporates re- 
cent experimental results on the compactification of porous aggre- 
gates, suggests that coagulative growth beyond millimetre sizes is 
difficult to achieve because compacted aggregates tend to bounce 
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rather than stick. These latter results raise serious questions about 
the validity of the incremental picture of planetesimal formation. 

A number of the difficulties faced by the incremental model 
have been known for many years, leading to searches for alterna- 
tive planetesimal formation scenarios. Gravitational instability of 
the dense dust layer that forms thro ugh vertical settling was sug- 
gested by iGoldreich & Ward! (Il973h as a means of forming km- 
sized planetesimals. Although this pathway is not now widely ac- 
cepted, recent models of planetesimal formation that involve col- 
lective effects, including self-gravity of regions o f enhan ced solids 
concentration, have been proposed. ICuzzi et al] ( 1200 Sl) have de- 
veloped a model in which mm-sized chondrules are concentrated 
in turbulent eddies, and then contract on the settling time of the 
solid grains under the action of self-gravity to form large 1 00 km- 
size d planetesimals. The comb ined action of the streaming instabil- 
ity jYoudin & Goodrriaij|2005l) . and trapping of metre-sized boul- 
ders in short-lived turbulent eddies, has been suggested as a means 
of in ducing planetesimal formation via gravitational collapse by 
Ijohan sen et al. ( 2007), leading to the direct formation of planetesi- 
mals that are as large as the asteroid Ceres {R^ ^ 500 km). Although 
these methods of forming planetesimals avoid the problem of grow- 
ing beyond the metre-sized barrier, they each have their own diffi- 
culties, and as yet no clear consensus of how planetesimals form 
has emerged. 

Once planetesimals have formed, runaway growth may lead 
to the formation of ~ 1000 km oligarchs on time scale s of ~ 10' yr 
jWetherill & Stewart|[T993l : iKenvon & Bromlevll2009l) . which then 
accumulate smaller bodies to become planetary embryos and cores 
durin g the oligarchic phase dlda & Makinol I1993I; iKokubo & Idd 
1 19981) . Rapid runaway growth of km-sized planetesimals requires 
that their velocity dispersion remains significantly smaller than the 
escape velocity of the accreting bodies, which for 10 km-sized ob- 
jects is ~ 10ms"'. A self-stining swarm of planetesimals embed- 
ded in a laminar disc can maintain a low velocity dispersion via 
gas drag and dynamical friction, but the situation in a turbulent 
disc may be very different due to stochastic forcing of planetesimal 
random motions by the turbulence. Planetesimal accretion that oc- 
curs at a rate determined by the geometric cross section may lead 
to planetary formati on time scales in e xcess of observed disc life 
times of a few Myr ( lHaischetalj200lh . 

A further constraint is that planetesimal collision velocities 
must remain below the catastrophic disruption threshold for col- 
lisional growth to occur. For planetesimals between the sizes 1 
- 10 km, disruption thresholds lie approximately in the interval 
10 ms"' < v„it < 60ms"', depending on the material strength 
and whether the planetesi mals are monolithic bodie s or rubble piles 
teenz & Asptiaug. 1999;,Stewart & Leinhardtl2009h . Clearly an ex- 
ternal source of stirring, such as disc turbulence, may not only pre- 
vent runaway growth, but may lead to the collisional destruction of 
planetesimals. 

Turbulence is generally believed to provide the anomalous 
source of viscosity required to e xplain the typical ~ lO'^^M ^yr"' 
accretion rates of T Tauri stars dSicilia-Aguilar et al.ll2004h . The 
most promising source of this is MHD turbulence driven by the 
magnetorotational instability (MRI, Balbus & Hawlev I991I). Both 
local shearing box and global simulations indicate that the non lin- 
ear, saturated state of the MRI is vigorous turbulence that is able 
to transport angular momentum at a rate that can exp l ain obser- 
vations of T T auri accret ion rates (j Hawlev et al." 1 9951 : lArmitag3 
ll998l : lHawle"vl [2001: Pa paloizou & Nelson 2003). But the MRI re- 
quires that the gas i s sufficiently ionised i n order for the linear in- 
stability to operate ( lBlaes&Balbuslll994h , and for non linear tur- 



bulence to be sustained jFleming et alj2000l) . In the absence of suf- 
ficient free electrons in the gas phase, the flow returns to a laminar 
state. 

Protoplanetary discs surrounding T T auri stars are cold and 
dense, leading to low levels of ionisation l lUmebavashilll983h . In 
the main body of these discs the prima ry sources of ionisatio n are 
expected to be external: stella r X-rays ([Glassgold et al.lll997h and 
possibly galactic cosmic rays JUmebavashi & Nakanol 198 1[). each 
of wh ich have a limited penetration depth into the disc. iGammid 
( 1 19961) presented a model in which the disc surface layers are 
ionised, and sustain a turbulent accretion flow, whereas the disc in- 
terior is a 'dead zone' where the flow remains laminar and minimal 
accretion takes place. This layered accretion disc model has been 
the subject of numerous investig ations over recent years which have 



examin ed the role of dust grains jSano et al .l2000t Ilgner & NelsoiJ 
l2006ah . gas-phase heavy m etals 1 Fromang et al.l 2002h . different 
chemical reaction networks (Ilgner & NelsonI 2006a ). the role of 



turbulent mixing dllgner & NelsonI l2006bl) ." and the Hall effect 
(Wardle 199^. 

MHD simulations have also been used to examine the dynam- 
ics of magnetised discs in which a finite el ectrical conductivity 
plays an important role. [Fleming et all ( |2000|) examined the satu- 
rated le vel of turbulence as a fu nction of the magnetic Reynolds 
number. iFleming & Stone! ( |2003|) examined disc models in which 
magnetic resistivity vari ed with height, s imulated the layered accre- 
tion model proposed bv lGammidd 19961) , and demonstrated that the 
dead zone maintains a modest Reynolds stress due to waves being 
excited by the active layers. More recent studies have coupled time 
dependent chemical networks to the MHD evolution using a multi- 
fluid approach, and have demonstrated that turbulent diffusion of 
charge carriers int o the dead zone can enliven it, in the absen ce of 
small dust grains dTumer et al]|2007l: lllgner & Nelsonll2008b . The 
presence of a significant population of small grains, however, re- 
duces the time for removal o f free charges from th e gas phase, and 
the dead zone is maintained dTumer & Sanoll2008h . 

Low mass planets and planetesimals embedded in turbulent 
protoplanetary discs are subject to st ochastic gravitational forces 
due to turbul ent density fluctuations dNelson & Papaloizoull2004l : 
lNelsonl200^ . The stochastic forcing leads to diffusion of the semi- 
maj or axes and excitat i on of t he eccentricity of the embedded bod- 
ies. iNelson & Gressell d2OI0l) - hereafter paper I - examined the 
dynamics of planetesimals embedded in fully turbulent cylindrical 
disc models without dead zones, and demonstrated that the exci- 
tation of the velocity dispersion will lead to collisional disruption 
of planetesimals of size < 10 km. They also demonstrated that 
very good agreement can be obtained between local shearing box 
simulations and global simulatio ns when shearing boxes of suffi- 
cient size are employed (see also lYang et al.ll200^ . In this paper, 
we extend the work presented in paper I and examine the dynam- 
ics of planetesimals embedded in vertically stratified discs, with 
and without dead zones, using shearing box simulations. A key re- 
sult is that the turbulent stirring of planetesimals is significantly 
reduced in discs with dead zones, possibly allowing for their con- 
tinued growth rather than destruction or erosion during collisions. 
As such, we propose that dead zones may provide safe havens for 
planetesimals, and are a required ingredient for the formation of 
planetary systems. 

This paper is organised as follows: We describe the physical 
model and numerical methods in Section |2] The presentation of 
the results is divided into four parts: In Section|3] we describe the 
general morphology and dynamics of the emerging hydromagnetic 
turbulence. Section |4] discusses how the resulting stochastic grav- 
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itational torques experienced by bodies embedded in the disc can 
be classified by their distribution and auto-correlation. Finally, in 
Sections[5]and[6l we analyse the temporal evolution in eccentricity 
and radial diffusion of a swarm of particles immersed in the flow. 
We discuss the implications of our results for planetary formation 
in Section[71 and we draw our conclusions in Section|8] 
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2 METHODS 

For the models reported-on in this pa per, we make use o f the 
second-order Godunov code nirvana-hi l lZieglerll200l |2008|) and 
solve the standard equations of resistive magneto-hydrodynamics 
(MHD) in a local, corotating Cartesian frame (x, y, z), which read: 



<9,p + V-(pv) 
5,(pv) + V-|pvv + p*- BB] 
a,B- Vx(vxB-?7VxB) 
V B 



0, 

IpQ. (^flxic - zxv) - pVO(z)z 

0, 

0, 



tQ} V, the total pres- 



with a static gravitational potential 0(z) 
sure p* = p + ^B-^, and where all other symbols have their usual 
meanings. The first source term in the momentum equation in- 
cludes the Coriolis ac celeration and tidal force in the shearing 
box a pproximation (see lGressel & Ziegleil20"07l : IStone & GardineJ 
I2OI0I for implementation details). For the shear rate, q, we apply 
the Keplerian value of -3/2. 



2.1 Model geometry and parameters 

The simulations are semi-global in the sense that we only con- 
sider a small horizontal patch of a protoplanetary accretion disc 
but include the full vertical structure. Our fiducial stratified simu- 
lation covers ±4 pressure scale heights, H, to account for hydro- 
static stratification under the assumption of an isothermal equation 
of state. To accommodate sufficiently wide MRI-active regions in 
the presence of a realistic dead zone, we use an increased box size 
of 5.5 scale heights on each side, i.e., a box with L. = 11 H. This 
provides an extra 3 to 4 above and below the nominal dead zone, 
whose vertical extent either side of the midplane is ~ 2H. 

Based on our previous work on permissible shearing box sizes 
for studying the dynamics of embedded bodies in turbulent discs in 
paper I, we chose a standard horizontal extent of 4 x 16 H. Due to 
the higher computational cost of the models including a dead zone, 
we adopt a somewhat smaller horizontal domain size of 3 x 12//, 
yet stil l large enough to allow for the e xcitation of spiral density 
waves jHeinemann & Papaloizou|[2009a| [bl). The box sizes and grid 
resolution for the different models are listed in Tab. \T\ 

Because the ionisation model used to compute the resistivity 
involves chemical rate equations, we need to specify a unit sys- 
tem to convert between computational units and physical units for 
the mass density and temperature. To allow for direct compari- 
son with other results, we scale our model such that it approxi- 
mates the wide ly used "minimum-mass" protosolar nebula model 
( lHavashilll98lh . More specifically, we place our local box at a ra- 
dius ro = 5 au and choose a disc aspect ratio h = H/r = 0.05, 
a column density of S = 135 gcm"^, a temperature T = 108 K, 
and an isothermal sound speed of c, = 667 ms"'. Using our def- 
inition of H, the equilibrium vertical density profile is given by 
piz) = po exp |^-z^/(2//^)j, where po is the midplane density. 



Table 1. Simulation parameters for the studied models. Columns 4-6 indi- 
cate the ionisation flux for X-rays (XR), short-lived radionuclides (SR), and 
cosmic rays (CR), respectively; open circles indicate nominal values, and 
a multiplication symbol followed by an integer gives the multiple of the 
nominal value adopted. 



2.2 Initial and boundary conditions 

For the initial magnetic field, we apply the standard zero-net-flux 
magnetic configuration B-{x) ~ sin(2;r x/L J, with an initial plasma 
parameter (ratio of thermal to magnetic pressure) ySp ^ 50. To avoid 
low Pp in the far halo region, we attenuate B^(z) with height, as to 
keep the plasma parameter roughly constant. At the same time, to 
obey the divergence constraint, we introduce a radial field /?j(x, z), 
which is subsequently sheared out into an azimuthal field as the 
system evolves. The described configuration has proven to produce 
a relatively smooth transition into developed turbulence, avoiding 
extreme field to pologies in the l i near g rowth phase. 

Following [Turner & San 3 ( l2008h . we furthermore impose a 
weak additional B- net field, with an associated midplane y6p = 
13900. Due to the assumed periodicity, the flux associated with this 
field is preserved for the duration of the simulations, and we found 
this to be a useful means of maintaining active MRI over long evo- 
lution times in the presence of a dead zone. Such a field may be a 
local residual of a large-scale field dragged in by the disc during its 
formation from the collapse of a magnetised molecular cloud. We 
also note that in the presence of a mean-field dynamo mechanism 
(as e.g. discussed inlOressel 201(1), one would expect a significant 
net vertical flux in any given local patch of the disc. 

For the velocity field, we use vertical boundary conditions 
which allow material to leave the box but prevent inflow. This has 
the advantage that the magnetic field can escape the box rather than 
pile-up nea r the domain boun dary, as is the case for a completely 
closed box ( IStoneetal.lll996l) . The magnetic field boundary condi- 
tions implement a rough approximation to an external vacuum field 
and are usually refeiTed to as "pseudo-vacuum" conditions. This 
means that the horizontal components are set to zero on the bound- 
ary, and only a vertical field component (with vanishing gradient) 
is allowed. We note that this type of boundary enforces a gradient 
in B, and By near the top and bottom of the domain, which results 
in an accelerated wind in this region, clearing material and em- 
bedded field structures from the magnetically dominated halo re- 
gion. To compensate for the mass loss associated with the outflow 
boundaries and to enforce a stationary background, we continually 
re-instate the initial density pr ofile in each grid cel l by means of an 
artificial mass source term (cf. iHanasz et al.l l2009h. 

In the simulations adopting a dead zone due to insufficient 
ionisation we prescribe a locally and temporally varying mag- 
netic diffusivity n. Som e e arlier simulations of dead zones by 
[Fleming & Stond l2003h and lOishi et"al] ilOOH) prescribed a static 
diffusivity profile derived from a physically motivated ionisation 
model. Here we further extend the realism of the diffusivity model 
by deriving a time-variable profile based on the gas column den- 
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sity, the detailed procedure of which is described in the following 
section. 

After our disc models have achieved a steady turbulent state, 
we introduce 25 planetesimals into the disc midplanes. We treat 
these planetesimals as non interacting test particles which evolve 
under the full 3D influence of the gravitational field of the turbulent 
disc. The disc gravity is softened on a length scale equal to the cell 
diagonal. 



2.3 The difTusivity model 

Co mparing different chemical reaction networks, ISano et al.lflOOOl) 
and lllgner & Nelson ' (2006a) found that the adsorption of free elec- 
trons and ions onto the surfaces of small dust grains plays an im- 
portant role in removing free charges from the gas phase. These 
authors find that even tiny mass fractions of micron-sized grains 
can significantly deplete charge carriers. 

Taking this dominance of small grains into account, we 
here refrain from following the detailed non-equilibrium chem- 
istry in combinatio n wi th a multi-fluid t r eatme nt as was done by 
iTumer et all ( l2007h and lllgner & NelsonI fcOOSh . Instead we adopt 
a simplified treatment of the gas-phase reactions. As a reason- 
able approximation to the more intricate modelling outlined above, 
we assume that recombination happens much faster than any dy- 
namical mixing timescale in o ur system, which seems warranted 
in the presence of small grains jSano et al ] |2000tlllgner&Nelsorj 
[200 6a b). As a first step towards more realistic description, we al- 
low the magnetic diffusivity, rj = //(x, t), to vary spatially in all 
three dimensions, and update it according to a look-up tabl e derive d 
from the reaction network in model4 o f lllgner & NelsonI ( l2006ah . 
In particular, we assume dust grains of size O.lyum and with den- 
sity 3 g cm"^ and use a dust-to-gas mass ratio of 10"^ . The assumed 
gas-phase Magnesium abundance is taken to be depleted by a factor 
10"'' compared to its solar value. The input parameters to the tabu- 
lated diffusivity are: (i) mass density, (ii) gas temperature, and (iii) 
the local ionisation rate All three quantities are evaluated on a 
per-grid-cell basis. To include the effects of external irradiation on 
(, we compute column densities to both the upper and lower disc 
surfaces. 



2.3.1 Ionisation sources 

There are several lik ely sources of ion i sing r adiation in the vicinity 
of a newly bom star. iTumer & Drakel (l2009h have recently studied 
how contributions from stellar X-rays, radionuclides, and energetic 
protons (interstellar cosmic rays, protons accelerated in the coronae 
of the star and disc) influence the shape and extent of a possible 
dead zone. We have implemented their ionisation model, including 
all ionisation sources. Because of uncertainties related to some of 
the proposed mechanisms, we take a conservative standpoint in this 
paper and focus on stellar X-ray irradiation (XR), and interstellar 
cosmic rays (CRs) as the prime sources of ionisation. 

Observations of young stars in Orion show that T Tauri star 
X-ray luminosities vary by approximately four ord ers of magni- 
tude, with a median value of Lx r - 2 x 10'" ergs"' l lGarmire et al.l 
I2OOOI) . Ilgea & Glassgoldl il99^ performed Monte-Carlo radiative 
transfer calculations of the ionisation rate in a standard protostel- 
lar disc model, adopting the above value for Lxr and assuming a 
plasma temperature of kT = 5 keV. Applying a fit to their results, 
we approximate the ionisation rate due to X-rays by 
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Figure 1. Vertical profiles for runs Dl and D2, in temis of the magnetic 
Reynolds number, Rm = O.H^ '7(z) ' , as given by our ionisation model. The 
thick line corresponds to the adopted initial model (with Z = 135 gcm~^). 
Thin lines on either side are for disc masses increased and reduced by fac- 
tors of two, and demonstrate the stron g dependence on the m ass column. A 
dead zone is expected for Rm S 10 (cf. lTumer & San 320081) . 



where is the position of our box in astronomical units, S„ and 
E/j, are the gas column densities above and below a given point, 
and ExR = 8.0 g cm"* is the absorption depth for the assumed spec- 
trum of X-rays. As detailed in [Umebavashi & Nakano (200i), we 
prescribe the following vertical attenuation of interstellar CRs illu- 
minating the disc surfaces: 



: 5x10" 



1 + 



A. 

2cR 



(2) 



Here, Sc r = 96 gem ^ is the cosmic r ay attenuation depth esti- 



mated by Umebavashi & Nakand (Il98lh . and the dots indicate the 
comsponding contribution from the second column density Sj. 

Because the numerically allowed diffusive time step decreases 
with the grid spacing squared, adequately resolved dynamical sim- 
ulations can become prohibitive in the presence of high diffusion 
coefficients. Given limited computational resources, we therefore 
decided to confine the dynamical span in the magnetic diffusivity 
77 to a reasonable range. We primarily do this by including an am- 
bient ionisation due to the decay of short-lived radionuclides (SR), 
which we further enhance by a factor of ten compared to the fidu- 



^xr ■ 



(1) 



cial value of <fsR = 3.7x 10" s" given in Turner & Drake C 
As can be seen in Fig.Q] where we plot the resulting 77 profiles, the 
ffoor in the magnetic Reynolds number due to the ambient level of 
ionisation is well below the expected threshold for the formation 
of a dead zone. This implies that we should not expect our results 
to be significantly affected by enhancing the effects of decaying 
radionuclides. 

We remark that, compared to previous work, our diffusivity 
pro files overall imply lo wer magnetic Reynolds numbers. Given 
that lFleming et al.l hoOOl) have found a critical Rm^ ^ lO'* for sus- 
tained MRI in a zero-net-flux configuration, we expect that the am- 
plitude of an external net-flux will have a significant impact on the 
level of turbulence. 



2.3.2 Sub-cycling 

Notwithstanding the moderate range in 77, the numerical constraint 
given by the diffusive propagation of information is one of the main 
limiting factors in our simulations. To avoid a degradation of the 
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numerical accuracy in the MRI-active region by an excessively low 
numerical time step in the dead zone, we choose to account for the 
two separate regimes in Rm. This is done by operator-splitting the 
diffusive term in the induction equation and applying a sub cycling- 
scheme for its updat^H- By doing so, we are able to integrate the 
non-diffusive part of the MHD equations with a longer time step, 
enhancing the accuracy of the solution in regions of low magnetic 
diffusivity. 



2.4 Code improvements 

Properly resolving the growing modes of the magneto-rotational 
instability (MRI) with Godunov-type codes has been found to de- 
pend on the reconstruction strategy used and on the ability of the 
Riem ann solver to capture the Alfvenic mode dSalsara & MeveJ 
l2oToh . To improve the representation of discontinuities in our finite 
volume scheme, we recently extended the nirvana-iii code with the 
Harten-Lax-van Leer Dis continuities (HLLD ) appro ximate Rie- 
mann solver introduced by iMivoshi & Kusand l l2005h . To guaran- 
tee the require d directional bias ing of the electromotive force in- 
terpolation (cf. lFlock et al.ll2O10l) . we h ave implemented and teste d 
the upwind reconstruction procedure of lGardiner & Stonel(l2005l) . 

To enhance the stability of our code in the strongly magne- 
tised corona, we gradually degrade the reconstruction order from 
second- to first-order accuracy near the vertical domain boundaries, 
thus avoiding undershoots in the hydrodynamic state variables in 
strong shocks. 



2.4.1 Artificial mass diffusion 

To fac ilitate the study of a low-beta disc corona, iMiller & Stond 
( l2000h have introduced the concept of a so-called Alfven speed lim- 
iter, circumventing prohibitively high signal speeds in low-density 
regions. Such a limiter can in principle be adopted for the approx- 
imate Riemann solvers we use. We chose, however, a different ap- 
proach and instead add an artificial mass diffusion term to the equa- 
tions of mass and momentum conservation. The diffusion coeffi- 
cient is chosen to lie well below the truncation error in the bulk 
of the domain. In grid cells where the density contrast exceeds a 
specified dynamic range, Cjy„, the coefficient is gradually adjusted 
according to 



1 + 1 10^^>» ^ 
Po 



-1/4 



(3) 



such that the grid Reynolds number approaches order unity - co- 
inciding with the stability limit for the explicit time integration of 
the diffusive term. The transition function is chosen in a way that 
the inverse operation can be efficiently implemented by means of 
consecutive square-root operations. We typically chose Cdyn = 5, 
resulting in four to five orders of magnitude in density contrast. 
We note that, because the diffusive fluxes are part of the finite vol- 
ume update, this approach does not violate mass conservation, and 
therefore avoids the problems of enforcing an artificial floor value 
in the density. In conclusion, our approach has proven to greatly 
benefit the overall robustness without noticeably affecting the solu- 
tion in the inter i or of the d omain. Similarly to the lim iters used by 
iMiller & Stond ( l2000h and ljohansen & LevinI ilOOA the Courant- 
Friedrichs-Lewy time-step constraint due to the fast magnetosonic 

' We typically cfiose a sub-cycling ratio of 5 - 8 to approximately match 
the restriction given by the high Alfven speeds in the halo. 
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Figure 2. Visualisation of the turbulent flow stracture. The colour coding 
indicates the vorticity, logiQ |V X u| of the perturbed flow u = V - qClxy. 
While the turbulent regions above and below the dead zone show strongly 
folded vortex structures, the flow pattern near the midplane is dominated by 
shearing waves. The contrast in the vorticity amplitude is ~ 100. 



mode is significantly alleviated, allowing for economical use of 
computing resources. 



3 DISC MODEL PROPERTIES 

In the following, we will present results from three different simula- 
tions (see Tab.[T]for details). To be able to make a quantitative com- 
parison with respect to the effects of a dead zone, we performed a 
fiducial stratified model with fully active MRI, referred to as model 
Al. In general, this model agrees well with the unstratified sim- 
ulations presented in paper I. For our standard d ead zone model, 
Dl, w e choose the reference ionisation rates of iTumer & Drak3 
l l2009h along with a dust-to-gas mass ratio of 1 : 1000, i.e., account- 
ing for a modest depletion of the smallest grains by coagulation and 
sedimentation. The resulting diffusivity profile is plotted in the left 
panel of Figure [T] 

The resulting dead zone in this model covers roughly ±2H, 
making it a reasonable proxy for what a realistic dead zone at 
5 au might look like. Gaseous nebulae around newly forming stars, 
however, are subject to substantial variations in X-ray irradiation 
jGarmire et al.ll2000h , so for our second dead zone model, D2, we 
increase the stellar flux by a factor of twenty. We find the average 
total thickness of the dead zone to be reduced by roughly 1.5 pres- 
sure scale heights in this case (see Fig.[T] right panel). 



3.1 Hydromagnetic turbulence 

One im portant result of the earl y simulations of layered accretion 
discs bv lFleming & Stond ( l2003h was that the dead zone, despite its 
name, retains a non-negligible level of Reynolds stresses, namely in 



© 2002 RAS, MNRAS OOO.fnfTS] 



6 Gressel, Nelson & Turner 





orbits 


(a) 


r 2 -2i 

cTr I cm s J 


Tf [2nQ. '] 


C^(Ax) [//] 


Cir(e) [Hlr] 


Al 


217 


0.0105 


0.45x10** 


0.30 


5.21x10^3 


2.68x10-3 


Dl 


224 


0.0038 


0.06x10** 


0.29 


4.72x10-'' 


1.54x10-* 


D2 


223 


0.0051 


0.13x10** 


0.27 


7.25x10-'' 


2.50x10'" 


Bl' 


505 


0.05 


0.74x10** 


0.32 


7.70x10-3 


2.77x10-3 



'torque corrected for 2D/3D evaluation (cf. Fig. 8 in iNelson & Gressel2O10l) . 
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Figure 3. Comparison between the box averaged turbulent Maxwell- and 
Reynolds stresses. Values are normalised by the gas pressure in the disc 
midplane, po- The time-averaged sums of these, i.e., the a parameters, are 
listed in column three of Tab.|2] 



the form of waves that are excited in the active layers. These resid- 
ual motions can clearly be seen in Fig.|2l where we have visualised 
the flow structure in terms of the logarithmic vorticity. The colour 
coding exhibits very distinct patterns in the two regions: while the 
MRI active layers show folded vortex features characteristic of de- 
veloped turbulence, the dead zone region is clearly non-turbulent, 
and dominated by sheared density waves with kn » k^. At the 
interface between the two zones, one can furthermore see struc- 
tures indicating some level of turbulent mixing into the dead zone. 
Because our model assumes that recombination on the surface of 
small dust grains is efficient (leading to a recombination time scale 
much smaller than the dynamical mixing time scale) this does not, 
however, have an effect on the extent of the dead zone. 

Our primary aim in this paper is to examine statistically the 
dynamical evolution of planetesimals embedded in turbulent discs 
with and without dead zones. In order to do this we require a quasi- 
stationary model of turbulence. As can be seen in Fig.|3] where we 
plot the time variation of the box-averaged transport coefficients, 
our models fulfil this requirement of a stationary mean with (ad- 
mittedly strong) supeiposed fluctuations. Because the bulk of the 
transport is associated with the MRI-active regions, the turbulent 
stresses only seem to depend weakly on the actual width of the 
dead zone. We remark that due to the different field configuration, 
the overall strength of the turbulence is somewhat weaker than in 
our previous non-stratified models presented in paper I. This is re- 
lated to buoyant loss of magnetic flux from the disc in stratified 
sim ulations, and is c onsistent with recent results in the literature 
(cf lDavis et"ai]|201(l and references therein). 



3.2 Spatio-temporal disc evolution 

Figure|4]shows space-time plots of the mean toroidal magnetic field 
By{z, t), as well as the total Maxwell stress for the fully active model 
Al. Probably the most striking features in this plot are the peri- 



Table 2. Overview of simulation results. Gravita- 
tional torques T represent the standard deviation o" 
of a normal distiibution (see Fig.|8], Coherence times 
Tc are from a fit to the torque ACF (see Fig. |9}. The 
random walk amplitudes Co-(A.y) and Cn-(e) refer to 
the dispersion in the displacement and eccentricity of 
a swarm of particles. 



odic "butterfly" patterns in the large-scale magnetic field, which 
are commonly observe d in this type of simulation (e.g. lOavis et al.l 
l2010l : lFlaig et alJboiOh . Because the upwelling field structures are 
not associated with a bulk motion of the flow, the likely explana- 
tion for the emer ging patterns is a mean-field dynamo as recently 
re-investigated bv lGressell | |2010|) . Even though we see quite strong 
fluctuations in the total Maxwell stress, the overall state is quasi 
stationary. Note that the two most prominent peaks of the total 
Maxwell stress (i.e. around 140-145 and 155-160 orbits) are (a) 
due to a strong mean field, and (b) confined to one half of the box. 
The situation of magnetic activity being dominan t in one half of 
the di sc has been observed before (see e.g. fig. 7 in lMiller & Stond 
I2000I) . and might be explained by equally strong quadmpolar and 
dipolar contributions. These cancel on one side and, at the same 
time, add-up on the opposite side, resulting in the observed lop- 
sided appearance. In fact, we cannot define a clear parity of the 
field in Fig. |4] which implies that the dipolar and quadmpolar dy- 
namo modes seem to possess similar growth rates. 

In the following discussion, we will focus on the vertical struc- 
ture of the disc and how it is shaped by the presence of a dead zone. 
For this we will show the time-averaget|l vertical profiles of vari- 
ous quantities. To demonstrate that this is a meaningful procedure, 
we provide space-time diagrams of three representative quantities 
(Fig. [5) for model Dl. A gain, the mean toroidal field generally is 
of irregular parity, i.e., neither quadmpolar nor dipolar symmetry is 
seen to prevail. Whether the concept of a global parity is relevant 
in the presence of a highly diffusive "insulating" layer remains a 
matter of discussion. It is interesting to note, however, that MRI 
channel modes - which are a non-linear solution to local net-flux 
MRI simula tions - represent g lobal modes even when stratification 
is included l lLatter et ani201(ll) . In fact, we see some indication of 
channels in our simulations, most prominently in the density (com- 
pare the lower panel in Fig. l5]with lLatter et al.ll2010l fig. 3). 

Quite remarkably, the number of field reversals related to the 
dynamo-cycle seems to be largely unaffected by the presence of the 
dead zone. As mentioned above, one main motivation in explaining 
the rising pattems by a mean-field dynamo, was the absence of a 
bulk motion near the midplane - implying t hat the pattern speed 
cannot be explained by the Parker instability jShi et al ]|20IOl) . The 
key in explaining the upward motion was a negative buoyant a ef - 
fect near the midplane ( Rudiger & Pipiij2000l ; Brandenburg|2008h . 
With the pattern in the disc halo being unchanged in the pres- 
ence of a dead zone, this picture possibly has to be a mended. 
Naive ly, for a positive a effect in the halo (as found bv iGressell 
bold) , one would not expect an outward butterfly diagram from 
classical theory - demo nstrating the need for non-linear f eedback 
both in the diagnostics feheinhardt & Brandenburg"2010') as well 
as in the modelling |Brandenbu rg, Candelaresi & Chatter] ee 200i; 
ICourvoisier, Hughes & ProctoilboiOl) . 

We note however, that because of the presence of the weak net 



If not stated otherwise, we average over the interval [20,210] orbits. 
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Figure 4. Spatio-temporal evolution of (a) mean toroidal field, (b) total Maxwell stress, for the fully active model Al. 
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Figure 5. Same as Fig.|4] but for the dead zone model Dl, and including an extra panel showing the rela tive density fluctuations (logarithmic). The horizontal 
field exhibits the same dynamo cycles in the MRI active layers. As already seen bv lTumer & SanolfcOOSb . the toroidal field leaks into the diffusively dominated 
midplane region, contributing to the stresses. Despite the absence of strong turbulent fields, density fluctuations reach a level of 10-20% within the dead zone. 



vertical flux in our current simulations, they are not strictly compa- 
rable to our previous work. Due to the net flux, a possible explana- 
tion of the unchanged migration direction might be the presence of 
prevailing channel modes. Th ese were show n to produce a strong 
negative helicity (cf. fig. 4 in lGresse il l201(3l) . compatible with the 
upward motion of the wave pattern^ 

Now focusing our attention to the actual dead zone region, we 
see that a significant amount of azimuthal (and in fact, radial) field 
diffiises in to the low-Rm region around the midplane . This leak- 
age has first been seen in simulations bv lTumeretallteOOTh and 



This was already suggested as an alternative scenario in the mentioned 
paper. 



iTumer & Sanol ( |2008b , which are very similar to ours and also con- 
tai n a weak net vertical field. In contrast to this, the simula tions 
bvlFleming & Stonelll2003h , IOishi. Mac Low & Menoul | |2007|) and 
lllgner & NelsonI ( |2008|) ." which apply a zero-net-flux configura- 
tion, do not show such an eff'ect. It seems peculiar, however, that 
the presence of a net flux should play a role in this context. As 
lllgner & NelsonI ll2008h illustrate in their fig. 2, the assumed diffu- 
sivity profiles o f iFleming & Stond ( l2003l) result in moderately high 
magnetic Reyn olds numbers near the midplan e . On the contrary, 
in the model of iTumer. Sano & DziourkevitchI ( |2007|) and also in 
our model (cf. Fig[T) the magnetic Reynolds number reaches much 
lower values of Rm < 10 for \z\ < H, leading to a much reduced 
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height z [H] 

Figure 6. Maxwell stresses for the different models. Thick lines indicate 
the total stress oc - (BrB^), thin lines represent the turbulent contribution 
cx - {B'^B'^), where B' =B - (B). The vertically integrated mass accretion 
rates (see rhs axis) for the models D2, Dl and Al, are 9.85 X 10-'*,7.82x 
10"**, and 1.36 X 10"^ M© yr"', respectively. Note the unexpected trend in 
the depth between models Dl and D2, as opposed to Fig.]?] 

diffusion time scale, possibly facilitating the leakage of large-scale 
field from the disc halo into the dead zone. 

Finally, we note that relative density fluctuations remain at a 
level of up to ten per cent within the dead zone - as seen in the low- 
ermost panel of Fig.|5] Episodes of lower perturbance (as indicated 
by lighter colours) are clearly correlated with strong azimuthal 
field, and notably with a negative Maxwell stress at the interface 
between the dead and active zones. As already mentioned above, 
we see indications of persistent channel modes (dark streaks above 
|z| ^ 2H), which are not unexpected given the magnetic configura- 
tion. It might be argued that such features are overly pronounced 
at the cumnt numer ical resolution, and migh t be destroyed by 
paras itic instabilities jPessah & Goodmanll2009l : lLatter et al .120091 
when the resolution is further increased. 

3.3 Disc vertical structure 

Figure |6] shows time-averaged vertical profiles of the R(f> compo- 
nents of the Maxwell stress tensor. Because we use an isothermal 
equation of state, resulting in a constant sound speed, Cj, these pro- 
files translate directly into a fractional mass accretion rate as in- 
dicated by the right-hand axis. The values on this axis are given 
by the expression HdM/dz = 3n H- p(z)T,,i, / p{z), which can be 
derived from the usual approximation for the mass accretion rate 
in a steady disc M = 3nvl,, but relaxing the assumption that the 
effective viscous stresses are independent of z- 

The vertically integrated mass flow rates generated by the 
magnetic stresses are 9.85 x 10"^ and 7.82 x 10"* Mq yr"' for runs 
D2 and Dl, respectively. This similarity is expected because the 
heights of the dead and active zones about the midplane in each 
model are similar (~ 1.5// and 2H, respectively, for the dead zones 
in models D2 and Dl). The mass flow rate in the fully active run A 1 
is 1.36x 10"' Mq yr"', nearly twice that observed in the model with 
the largest dead zone. Given that the bulk of the mass in the disc 
lies within the dead zone, where the Maxwell stresses for models 
D2 and Dl are factors of ~ 10"^ and 10"^ smaller than for model 
Al, the bulk of the transport in models Dl and D2 happens in the 
active layers because the stresses (and effective viscosity coefficient 
a) remain large there. 

Figure |6] shows that the Maxwell stresses in the midplane 
regions are dominated by large-scale fields in runs Dl and D2 
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Figure 7. Same as Fig. |6] but for the Reynolds stress (pvuSv^). The 
vertically-integrated mass accretion rates are 2.08 X 10"*, 1.39 X 10-**, and 
3.64 X IQ-* Mq yr"' for models D2, Dl, and Al, respectively. 



- even more so than in the halo regions of the fully active run 
A l. This result agrees b roadly with recent simulations (cf. fig. 4 
in iTumer & Sanoll2008h which included similar physical effects. 
While the turbulent contribution to the Maxwell stress (thin lines) 
falls-off to - lO"' in the dead zone for both models Dl and D2, 
one can clearly see that the effect of the greater dead zone width 
in model Dl is compensated for by a shallower valley in the total 
Maxwell stress (thick lines). The inverse top-hat shape of the pro- 
files is likely related to large-scale field diffusing into regions of 
intermediate Rm, as seen in the topmost panel of Figure |5] Given 
that the transition region is closer to the low-/? halo for model Dl 
(with a wider dead zone), it seems plausible that the strength of the 
regular field is larger in this case. Clearly, the contribution due to 
the fluctuating fields becomes negligible in both cases. 

Considering that, for model D 1 , the floor value of the Maxwell 
stress is of the same magnitude as the Reynolds stress (cf. Fig.|7]l, 
such large-scale stresses might play an important role in mediat- 
ing mass transport within the dead zone - making a more detailed 
future study of the discussed trend worthwhile, particularly in the 
context of global disc models, where it may be possible to observe 
large scale accretion flows in dead zones. 

In Figure [7] we plot vertical profiles of the time-averaged 
Reynolds stress. As expected, the lines agree very well in the MRI- 
active regions. Comparing the models Dl and D2, we now see a 
correspondence between the dead zone width and depth, as ex- 
pected. It appears that larger and more massive active zones lying 
above and below the dead zone are more able to excite higher am- 
plitude sound waves which are able to propagate into the dead zone, 
generating a correspondingly larger Reynolds stress there. 

Based on two runs with different (zero net flux) field topolo- 
gies, [Fleming & Stoii^ ( l2003l) concluded that "the Reynolds stress 
in the midplane does not seem to depend on the size of the dead 
zone but rather the amplitude of the turbulence in the active layers". 
This is clearly not the case for our models Dl and D2, which show 
quite different Reynolds stresses in the dead zone, while having 
almost identical levels of Maxwell stresses - both in the total and 
fluctuating part - within the active region (see Figs.|6]and|7] respec- 
tively). While we reckon that the larger Reynolds stress is related 
to the larger mass of the active layer in model D2, this discrepancy 
may also be due to the presence of large-scale fi elds diffusing into 
the mi dplane region in our simulations. Note that [Fleming & Stond 
( l2003h did not mention such fields in their discussion. 
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Figure 8. Torque histograms for models Al, Dl and D2. The grey hnes indicate fitted normal distributions, for which standard deviations are given (note the 
difterent scale on the abscissae). While the torques within the dead zone are well approximated by a Gaussian, the distribution in model Al develops some 
deviations for large torques. 



4 CHARACTERISING THE GRAVITATIONAL FORCING 

As first suggested by the global simulations of lNelsoiil ( l2005h . the 
stochastic gravitational forces associated with density fluctuations 
from developed MRI turbulence have the potential to limit severely 
the growth of planetesimals. In paper I, we have confirmed this 
original finding in the framework of global simulations and local 
simulations with large enough box sizes, but in which the vertical 
component of gravity was neglected. In the following section, we 
extend this work to the case of stratified discs with and without 
a dead zone. To facilitate a direct comparison, we here focus on 
the case of a moderately strong net vertical field. In this respect, 
the results should be interpreted as upper limits. To obtain lower 
limits for the torque amplitude, models with nominal dead zones 
and applying a zero-net-flux configuration will be required. 

4.1 Stochastic torques 

Figure[8]shows distribution functions of gravitational torques along 
the y-direction in our simulation box. The histograms are computed 
from t he time series of a set of 25 particles in each simulation run. 
Unlike lOishi et alj ( l2007h , we do not observe transients in the time 
series, and the torques are indeed quasi-stationary in the interval 
[20,210] orbits, which has been used to obtain each histogram. As 
required for a simplified treatment in terms of a stochastic Monte 
Carlo model, the histograms are well represented by a normal dis- 
tribution. Notably, there are moderate power-law tails in the fully 
active run Al, indicating some level of non-Gaussian fluctuations 
(also cf. fig. 7 in paper I). The amplitude of the torques in the 
difl"erent runs roughly scales with the square-root of the midplane 
Reynolds stress. This arises because the density fluctuations scale 
linearly with the velocity fluctuations, and coincid es with the scal- 
ing found in section 3.3 of lBaruteau & LinI ( I2OI0I) . 

Apart from the torque amplitude, the ability to stochastically 
amplify particle motions depends on the degree of coherence in 
the fluid patterns. Building on our work from paper I, we here ap- 
ply the same formalism and study temporal torque autocorrelation 
functions (ACFs). We use the fitting formula introduced in section 
3.3 of paper I, i.e. 

5 r(T) = [{l-a) + a cos(2;r w t)] e-^'^' , (4) 

with free parameters, a, a>, and t^. As discussed in paper I, we as- 
sume two components to better fit the shape of the ACF: (i) an ex- 
ponential decay representing the truly stochastic part of the torque 
time series; (ii) a modulation due to "wavelike" behaviour, i.e., a 



model Al model Dl 
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Figure 9. Torque autocorrelation function (ACF) for models Al and Dl. 
The black lines indicate model fits (see text), showing that motions are co- 
herent over approximately one third of the orbital time scale. The ACF for 
model D2 is very similar to Dl and is hence not shown. 



negative dip in the ACF reflecting the fact that density enhance- 
ments swashing over the region of influence will produce torques 
of opposite signs when approaching and receding the test particle. 
Such a feature is also seen in the c orresp onding Figs. 9 and 15 of 
lOishi et al] 1*2007') and I Yang et al.l ( |2009|) . respectively, complicat- 
ing the interpretation of the results. 

In paper I, we have discovered certain trends in the ampli- 
tude of the latter effect, which are probably related to the finite size 
of the computational domain. Because the periodic box in some 
respect acts as a resonator, the wavelike component probably ap- 
pears too pronounced in local simulations. This demonstrates that 
one does well keeping in mind the li mitations of local boxes i n 
terms of the allowed dynamics (also see lRegev & Umurhanll2008l) . 
The parameter of interest, of course, is the correlation time of 
the stochastic perturbation. While waves represent ordered motion, 
it is this stochastic part that ultimately produces the random-walk 
amplification of the particle motions. 

In the following, we rely on the approach taken in Eqn. l|4j 
providing an effective deconvolution of the two effects. The result- 
ing fits for models Al and Dl are shown in Fig. |9] Unlike model 
B I (cf. fig. 9 in paper I), model Al shows a moderate level of mod- 
ulation. Given the horizontal box dimensions of 4x 16H, this is 
somewhat unexpected as we observed a consistent trend towards 
weaker modulation for larger box sizes in paper I. This being said, 
the fitted coherence time of T(. = 0.30 in model Al is in striking 
agreement with our earlier unstratified model. The torque ACF of 
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Figure 10. Comparison of 
three different model ACFs as 
given by Eqn. jj)- All curves 
have the same coherence time 
Tf = 0.3 and yet show quite 
dift'erent zero-crossings. 



the dead zone model D 1 shows a comparable amplitude in the mod- 
ulation and has an almost identical value of = 0.29. 

With respect to t he artificial forcing func tion used in the 2D 
global simulations of ISaruteau & LinI l l201(]t) . we note that our 
value for = 0.3 coincides with the one stated for their "reduced 
lifetime" case. Comparing our ACF with the one in their Fig. 2, 
there appears, however, to be a discrepancy of a factor of two in 
the first zero crossing. This can possibly be explained by a different 
modulation of the exponential decay as exemplified in Figure [TO] 
All curves have a common = 0.3. Moreover, the first two curves 
have a 60% wavelike modulations with frequency u) as indicated, 
while the third is a simple exponential. Owing to the confusion with 
respect to the definition of the coherence time, it seems worthwhile 
to further study how the spectrum of modes affects the apparent 
coherence. Neverthele ss, given the actual type of forcing used in 
lBaruteau&Linl ( l2010l) , i.e. a spectrum of wavelike motions, one 
may be surprised how similar the ACFs indeed look. Going back to 
Fig. (2] it however becomes obvious that it is not the turbulence we 
need to approximate but only the spiral density waves with uncor- 
related phases that it creates. 

The coiTesponding ACF for model D2 is omitted here, as it 
is very similar to Dl. The fitted coherence time for model D2 is 
Tc = 0.27, in good agreement with the other two runs. We therefore 
conclude that the presence (and depth) of a dead zone does have 
little influence on the temporal torque statistics and merely affects 
the amplitude of the stirring process. 



5 ECCENTRICITY STIRRING AND PLANETESIMAL 
ACCRETION 

To study the excitation of the velocity dispersion and the radial 
diffusion of embedded boulders, planetesimals, and protoplanets, 
in each run we integrated the trajectories of a swarm of 25 non- 
interacting test particles subject to perturbations from the gravita- 
tional potential of the gas disc. Unlike in paper I, we here focus 
on larger bodies and neglect the effects of aerodynamic interaction. 
This is warranted for solids with radii above ^ 100 m (cf. paper I 
where the results for planetesimals with sizes in the range 1 00 m - 
10 km were found to be similar for run times on the order of a few 
hundred planetesimal orbits). The upper limit for the size range 
considered here is given by the constraint that perturbations of the 
disc remain weak such that spiral wave excitation and gap-opening 
can be ignored. This is the case for planetesimals and small proto- 
planets, where feedback onto the disc can be ignored. 

Although we do not consider in detail the dynamics of metre- 
sized boulders that are strongly coupled to the gas via drag forces in 
this paper, we recall from paper 1 that such bodies quickly achieve 
random velocities very similar to the turbulent velocity field of the 
gas. The midplane r.m.s. turbulent velocity for the fully turbulent 
model Al is found to be vvi, = (76.6 ± 6.5) ms"', and for model 
Dl vwb = (17.8±2.6)ms"'.FormodelD2 vturb = (28.5±3.6)ms"'. 
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Figure 11. Random-walk eccentricity growth for the dispersion, (T(e), of 
sets of 25 particles in each run. The expected V' behaviour is indicated 
by the dashed line. The rhs axis gives the coiTesponding radial velocity 
dispersion, leading to disruptive relative velocities for the fully active run 
Al. The amplitudes for the dead zone runs are lower by a factor of 10-20. 



We now consider the dynamics of larger planetesimals. The 
excitation of the eccentricity (and equivalently the radial velocity 
dispersion) can be seen in Fig. [TT] where we plot the r.m.s. val- 
ues of these quantities averaged over the planetesimal swarms as a 
function of time for the runs Al, Dl and D2. The time history is 
consistent with a random-walk behaviour, indicated by the dashed 
line representing an e{t) ~ -sft dependence. Expressing the time 
evolution of the eccentricity according to 



e{t)lh = CAe)^U 



(5) 



where h = ///r, and the time is measured in local orbits, we obtain 
the following values for the amplitude factors, Co-(e) (which are 
also fisted in Tab.Hi: 2.68 x W-\ 2.50 x 10-^ and 1.54 x 10"^ for 
models Al, D2, and Dl respectively. Despite the different levels 
of turbulence (as reflected in both a and F), the values in Tab. |2] 
show that the stirring amplitude of model Al agrees quite well with 
the unstratified model from paper I. The amplitudes from the dead 
zone runs D2 and Dl are lower by a factor of ^ 11 and ^ 17, 
respectively, showing that the presence of the dead zone has a clear 
and measurable effect on the strength of eccentricity stirring. 



5.1 Long-term evolution 

As we have only been able to run our simulations for about 200 
orbits, it is important to consider the long term evolution by exam- 
ining the magnitude of the equilibrium eccentricity that is expected 
to arise through the balance between turbulent eccentricity excita- 
tion and gas drag and/or collisional damping (Ida et al. 200_^. The 
question of whether or not collisional growth of planetesimals may 
occur in the different disc models is determined by the velocity dis- 
persions obtained relative to the disruption/erosion velocity thresh- 
olds, and this point is addressed in the discussion that follows. 



5.1.1 Gas drag damping 

We have not included gas drag or collisional damping in the sim- 
ulations presented in this paper, but gas drag was included in the 
runs presented in paper 1, and they showed that over run times of 
a few hundred orbits the gas drag has little influence on the eccen- 
tricity for planetesimals with radii R„ ^ 100 m. We adop t the usual 
expression for the gas drag force (IWeidenschiUing|ll977h 
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Table 3. Results for the equilibrium velocity dispersions obtained for each model as a function of the planetesimal sizes, and as a function of the damping 
mechanism. Also tabulated is the time required to grow to the smaller of the equilibrium velocity dispersion values given for each model each size. 



Fdmg = -ConRpP |Vp - Vgl (Vp - Vg) 



(6) 



where Cd is the drag coefficient. For making simple estimates 
we take the value Cd = 0.44, appropriate for larger planetesi- 
mals for which the Reynolds number of the flow around the body 
Ke > 800. The eccentricity growth over time is given by (|5) where 
the time, t, is measured in local orbits. Following the same proce- 
dure adopted in paper I, and calculating the equilibrium velocity 
dispersion by equating the eccentricity growth time scale and the 
gas drag-induced damping time scale, we obtain 



Vdis 



[ 10'CdP 



(7) 



where Rp and Pp are the planetesimal radius and density, and is 
the Keplerian velocity. Note that the form taken by l|7j now assumes 
that the system evolution time given in l|5j is expressed in seconds. 
The equilibrium velocity dispersions obtained using (|7) for each 
disc model and planetesimal size are listed in the third column of 
Tab.H 

5.1.2 Collisional damping 

The efl'ects of collisional damping, due to inelastic collisions be- 
tween planetesimals, may be estimated by noting that the colli- 
sional damping time scale is given by the product of the collision 
time scale and the number of collisions required to damp out the 
kinetic energy of random motion for a typical planete simal. This 
approach is similar to that adopted by llda et al.r ( l2008h . The colli- 
sion time (neglecting the effect of gravitational focusing) is given 
by 

1 



npnRjvi,sp ' 



(8) 



where ;ip is the number density of planetesimals. This may be ap- 
proximated by Hp = Zp/(2//p;77p), where Hp = vjisp/fij. and Zp is 
the mass surface density of planetesimals. The coefficient of resti- 
tution is given by Cr = I'f /vj, where Vj and I'f are the initial and final 
velocities associated with a collision. The number of collisions re- 
quired to damp I'disp is therefore 

1 



1-Cr 



(9) 



Combining Equations ^ and ^ we obtain the collision damping 
time: 

S^pPp / 1 



' c-daiTip 



1-Cr 



(10) 



Equating the collisional damping time with the eccentricity growth 
time associated with (O - and given explicitly by equation (16) in 
paper I - yields an expression for the equilibrium velocity disper- 



I'disi 



4RpPpCl(eWv] 
W-Lpflk 



1-Cr 



(11) 



The value to be adopted for the coefficient of restitution, Cr, 
depends on the material composition, the impact velocity, and 
whether or not the planetesimals may be considered to be mono- 
lithic bodies or rubble piles held together by self- gravity alone. 
Adopting the velocity dependent formula from [Bridges et al.l 
l ll984l) 



Cr = MIN 



0.34 (-^V 
\ 1 cm s ' / 



1 



(12) 



we find that Cr - for the values of the radial velocity dispersion, 
I'disp, that we obtain in the simulations. We therefore adopt Cr = 
when estimating the equilibrium i'disp. The equilibrium velocity dis- 
persions obtained using dllb for each disc model and planetesimal 
size are listed in the fourth column of Table |3] In obtaining these 
values we assume that Zp = 4x (S/250) - i.e. the surface density of 
solids is a factor of 1 /250 lower than the gas surface density, but is 
augmented by a factor of 4 beyond the ice-line. We further assume 
that all disc solids are incorporated within planetesimals of size ^p 
for each size that we consider. 



5.1.3 Disruption velocity thresholds 

The consequences of the equilibrium velocity dispersions obtained 
from equations ([7} and dllt for planetary growth can only be deter- 
mined by comparing them with the disruption/erosion thresholds 
for co lliding bodies jBenz & Asphaug|[l999l : IStewart & Leinhardl 
l2009l) . and with the esc ape velocities associated wit h the planetes- 
imals. In a recent study, IStewart & LeinhardJ (|2009|) present a uni- 
versal law for collision outcomes in the form 
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Mir ^ 

M,„, 2 2d ' 



(13) 



where M]^ is the mass of the largest post-coUision remnant, Mta, 
is the total mass of the colliding objects M[ + M2, and 2r is the 
reduced mass kinetic energy normalised by the total mass 



Qr 



2Mi 



-vt. 



(14) 



Here V\ is the impact velocity. For accretion to occur during a col- 
lision between equal sized bodies, we require M\JMi„t > 1 /2, or 
equivalently Qr/Qo < 1, such that 2d is the collisional disruption 
or erosion threshold. 

The value of 2d is sensitive to factors that influence the energy 
and momentum coupling bet ween colliding bodies (e.g. i mpact ve- 
locity, strength and porosity). IStewart & Leinhardtl(l2009h fit results 
from their numerical simulations and data in the literature using the 
expression 



(2-3,im) 



(15) 



2D = M:r'""*''+'?G^ir]^,' 

where q^, qo, <pm and are parameters. R12 is the spherical ra- 
dius of the combined mass, assuming Pp = 1 gcm"^. Since we use 
Pp = 2gcm"^, and consider collisions between equal-sized bodies 
only, we have Rn = The first term on the right hand side 

of J15t represents the strength regime, while the second term rep- 
resents the gravity regime. Small bodies are supported by material 
strength, which decreases as the planetesimal size increases. On 
the other hand, gravity increases in importance with growing plan- 
etesimal size - with the transition between the strength a nd gra vity 
regimes occurring at R^ x 100 m. IStewart & Leinhardti ( |2009|) de- 
rive the following values for the above parameters for weak aggre- 
gates (weak rock): = 500, qa = lO"'', /Jm = 0.4, 0,^ = 7. For 
stron g rocks they use basalt l aboratory data and modelling results 
from lSenz & Asphaud ( Il999h to obtain: q, = 7 x qo = 10 ^ 
yUm = 0.5, ipm = 8. In the limit of large planetesimals, the disruption 
curve can be best fit by their results for colliding rubble piles, for 
which 2d = 1.7 X 10"*i?p. Given the uncertainties associated with 
the structure and material strength of planetesimals, we tabulate the 
disruption velocities, Vcit, for both weak and strong planetesimals 
in the fifth and sixth columns of Tab.[3] For iJp ^ 10 km, we tabulate 
the disruption velocities for rubble piles in the seventh column. 



5.2 Model Al 

We now consider the results for the fully turbulent model Al. The 
equilibrium velocity dispersions coiTesponding to gas drag damp- 
ing are listed in the third column of Tab.|3]for 100m, 1 km, 10km 
and lOOkm-sized planetesimals (assuming a density pp = 2 g 
cm"^). These values of vjisp are smaller than those presented in pa- 
per I by approximately a factor of two, because in that paper we 
adopted a slightly more massive disc model, a larger planetesimal 
density pp = 3 g cm"^^, and we utilised cylindrical disc models that 
give rise to enhanced stirring of planetesimals relative to the full 
3D simulations considered here. 

Comparing the values of v^isp in the third and fourth columns, 
we see that for all planetesimal sizes damping due to gas drag dom- 
inates over collisional damping, so it is the values in the third col- 
umn that are closest to the true equilibrium values obtained when 
all sources of damping act simultaneously. We see that Vdisp > I'ciit 
for planetesimals composed of both 'weak' and 'strong' rock for 
Rf < 10 km, and for R^ = 100 km we see that Vdisp exceeds the 
disruption/erosion threshold for rubble piles. We conclude that if 



planetesimals reach their equilibrium values for vjisp, then mutual 
collisions will lead to their destruction. 

Estimated times for v^^sp to grow to the equilibrium values via 
a random- walk (vjisp ~ ^/t) are listed in the eighth column of Tab.[3] 
The most favourable models for the incremental formation of plan- 
etesimals at 5 au in a laminar disc suggest formation times of a 
few xlO^'yr JWeidenschillind [20001) . Formation times in a turbu- 
lent di sc exceed this beca use the vertical settling of solids is re- 
duced terauer et al1l2008h . The random velocity growth times in 
Tab. [3] thus indicate that if planetesimals were able to form incre- 
mentally in the disc that we simulate, then they would always have 
a velocity dispersion equal to the equilibrium value. But, the time 
to reach the disruption velocity v^^n = 15.2 ms"' for strong 1km 
aggregates is only 800 yr, much shorter than the formation time. It 
is clear that planetesimals cannot fomi and grow by a process of 
collisional agglomeration in a fully turbulent disc similar to model 
Al. 



5.3 Models Dl and D2 

We first discuss model Dl, which has the deeper dead zone of the 
two dead zone simulations. The values of v^jsp listed in Tab.[3]show 
that for planetesimals with Rj, < 10 km collisional damping domi- 
nates gas drag in setting the equilibrium velocity dispersion. Only 
for Rp = 100 km is gas drag more important. The transition from 
gas drag dominated damping in model A I to collisional damping in 
model Dl occurs because of the different functional dependencies 
on the eccentricity excitation coefficient, C,^(e) in equations Q and 

The equilibrium vjisp for all planetesimal sizes lie between the 
disruption/erosion thresholds, v^rit, for weak and strong aggregates. 
For the larger planetesimals with i?p = 10 and 100 km we see that 
^'disp < Vcrit for rubble piles. These results suggest that collisions 
between planetesimals in a dead zone can lead to growth rather than 
destruction, with the outcome depending on the material strength. 

The times required for vjisp to reach the equilibrium values are 
listed in the eighth column of Tab. [3] Given that the time required 
for incremental formation of km-sized planetesimals in a disc with 
modest turbulence is likely to exceed 10^ yr at 5 au, the time scales 
for the growth of v^isp suggest that planetesimals growing through 
a process of particle sticking will always have velocity dispersions 
close to the equilibrium values. The fact that these values are below 
the disruption thresholds for strong aggregates suggests that colli- 
sional growth may still be possible, albeit at a slower rate than in a 
laminar disc. 

If km-sized planetesimals can form, further accretion is 
normally expected to arise via runaway growth, leading to 
the formation of ~ 1000 km sized oligarchs within ~ 10^ yr 
jWetherill & Stewart! Il993l ; iKenvon & Bromle\i l2009h . Runaway 
growth requires the velocity dispersion to be lower than the es- 
cape velocity from the accreting bodies, and for an internal density 



2 g cm""* this is given by v^s^ = 10 (j^) ms"'. A typical 
10km planetesimal in model Dl forming via incremental growth 
will be excited to v^i^p > 10 ms"' during its formation, so turbu- 
lent stirring will prevent runaway growth from occuning for bodies 
of this size. Because I'disp remains below the disruption threshold 
for rubble piles, however, collisions can still lead to growth, but 
at a substantially slower rate than occurs during runaway growth. 
We discuss the implications of our results for the more rapid p lan- 
etesimal forma t ion sc enarios presented by ICuzzi et al. 1 (l2008l) and 
Ijohansen et alj ( l2007h in Section|7] 

Similar conclusions may be drawn from the results of model 
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Figure 12. Dispersion of the radial displacement Ax of a swarm of 25 test 
particles for run Al. The growing spread in the particles' separation is well 
approximated by a random-walk as illustrated by the over-plotted curve. 
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Figure 13. Model comparison for the random-walk behaviour shown in 
Fig. [12] above, but plotted logarithmically. The second axis gives the ab- 
solute dispersion Aa at r = 5 au. Evolution scaling with the square-root of 
time is indicated by the dashed line. 



D2 as we have drawn for model Dl. But, having a deeper dead 
zone, as in model Dl, is clearly preferable for planetesimal forma- 
tion due to the fact we find that the velocity dispersion is smaller in 
that case. 



6 RADIAL DIFFUSION OF PLANETESIMALS 

It was shown in paper I that the radial drift of planetesimals of size 
i?p > 10 km over typical disc life times due to gas drag is small. In- 
stead, the diffusion of planetesimal semimajor axes is dominated by 
the stochastic gravitational torques exerted on the planetesimals by 
the turbulent disc. This implies that material of common origin and 
composition (e.g. planetesimals which form at a particular radial 
location in the disc), in time will spread out over a certain region in 
semimajor axis. The degree of spreading is determined by the disc 
life time and the strength of the stochastic torques. 

In paper I, we examined the degree of spreading induced by 
the stochastic torques in a fully turbulent disc model similar in 
mass to the minimum mass solar nebula model, and showed that 
over a putative disc life time of 5 Myr, planetesimals embedded in 
the current asteroid-belt region would spread inward and outward 
over a distance of ~ 2.5 au. Such a conclusion sits uncomfortably 
with the relatively low vola tile content of the inner solar system 
planets l lO'Brien et ai]|2006h . which would presumably be higher 
if substantial volatile-rich material from beyond the snow-line had 
migrated inward during the planet-forming epoch. Additionally, al- 
though the picture of how different taxonomic classes of asteroids 
are distributed as a function of heliocentric distance h as become 
more co mplex in rece nt years ('Moth e-Diniz et al .l2003h . the earlier 
work of iGradie & Te desco ( 1982) shows a clear trend of increas- 
ing volatile content as a function of heliocentric distance. In this 
latter study, the distribution of taxonomic types is consistent with 
a picture of in situ formation, followed by radial diffusion over a 
distance of ^ 0.5 au. In the following discussion we judge that our 
models have had a measure of success if they do not violate the 
constraints implied by this picture. 

In Fig.[T2l we plot the evolution of the r.m.s. of the radial dis- 
placement of a swarm of test particles, and we find the resulting 
curve to be consistent with a random-walk. As mentioned above, 
we neglect dynamical interactions amongst the particles. This im- 
plies that self-stirring due to close encounters is not accounted for, 
and the given values have to be regarded as lower limits. 



In Sect.O we have demonstrated that the gravitational torques 
acting on particles can be represented by a normal distribution, and 
their temporal correlations possess a finite coherence time. These 
two properties allow us to estimate particle diff"usion based on a 
Fokker-Planck e quation. As dis cussed in detail in section §3.7 of 
paper I (see also iJohnson et al. I l200a) . the natural variable to de- 
scribe this process is the specific angular momentum j. The time 
scale for diffusion of particle angular momenta due to stochastic 
torques is 



(16) 



where Dj is the difl"usion coefficient and A j is the change in j. The 
diff'usion coefficient Dj may be approximated by Dj = crj.T^, where 
CTp is the r.m.s. of the fluctuating specific torques, and is the 
congelation time associated with the fluctuating torques. 

The change in specific angular momentum obtained after an 
evolution time of t is given by 



(17) 



Noting that small changes in specific angular momentum are re- 
lated to small changes in the semi-major axis, a, by the expression 



Aa 
2a 



the change in semimajor axis over an evolution time t is: 



Aa : 



./ 



(18) 



(19) 



We can now examine the degree of agreement between the 
level of particle diffusion observed directly in the simulations, and 
presented in Fig.[T3] and predictions based on ( 119) . The radial lo- 
cation of the shearing box in our simulations is assumed to be 5 au, 
and simulation run times are t =^ 200 local orbits. The value of Dj 
for each model may be computed using the expression Dj = ct^Tc, 
and values for <Tr, expressed in c.g.s. units, may be read off Fig. [8] 
for each model. The corresponding estimates of are listed in 
Tab.|2](along with the r.m.s. specific torques, cti-). 

We note that the time evolution of the r.m.s. radial displace- 
ment obtained in each model, cr^^., has been fitted using the expres- 
sion 



o-A,/// = a(Ax) V? 



(20) 



where H is the local disc thickness, the time t is measured in units 
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of local orbits, and the coefficients Ca-{^x) are tabulated in Tab.|2] 
After 200 orbits, the radial diffusion given by ( 120b corresponds to 
a change of semimajor axis Aa = 0.0184 au for model Al, which 
can be verified by inspection of Figure [13] The predicted change 
from M9\ is Aa = 0.0122 au, approximately 30% smaller than the 
observed value. 

For model Dl, we observe a change Aa = 0.00167 au, and 
predict a change from ( 119) of Aa = 0.00161 au, giving excellent 
agreement. For model D2, the observed change in semimajor axis 
is Aa = 0.00256 au, which is approximately 25% smaller than the 
predicted value of Aa = 0.00336 au. Taking these results overall, 
and noting that the use of only 25 particles in the simulations im- 
plies minimal sampling errors at the 20% level, we consider the 
agreement between the observed and predicted levels of diffusion 
to be very satisfactory. 

Examining the longer term evolution, we use ( 120) to calcu- 
late the level of diffusion expected over an assumed disc life time 
of 5Myr. For the fully active model, Al, we obtain Aa = 2.9 au, 
which is in good agreement with the result obtained in paper I. For 
model Dl, we obtain Aa = 0.26 au, and for model D2 we obtain 
Aa = 0.40 au. It is clear from these values that in a vertically strat- 
ified disc sustaining MHD turbulence without a dead zone, radial 
diffusion is predicted to be too large to be consistent with solar 
system constraints. For a disc with a significant dead zone whose 
height extends either side of the midplane by a distance ~ 2H, how- 
ever, we see that the degree of radial mixing is strongly diminished, 
and the results are consistent with the d istribution of asteroidal tax- 
onomic types dGradie & Tedescdl 19821) . 



7 DISCUSSION 

We now attempt to summarise our results and present a coherent 
picture of how the dynamics and growth of planetesimals are af- 
fected by turbulent stiixing in discs with and without dead zones. 
We frame our discussion in the context of the three planetesi- 
mal formation models discussed in the introduction: incremen- 
tal gro wth through particle sticking over time scales exceeding 
10^ yr ( IWeidenschilling|2000lJBrauer et al.l2008b : concentration of 
chondmles in turbulent eddies, followed by gravitational contrac- 
tion on the chondrule settling time - typically ~ 10' - 10^^ or- 
bits dCuzzi et"ai]|2008l) : gravitational collapse of dense clumps of 
metre-sized bodies formed in turbulent discs through a combina- 
tion o f trapping in local pre ssure maxima and the streaming insta- 
bility jjohansen et al.ll20()7n . Although there are significant prob- 
lems with the incremental growth picture, in this paper we are inter- 
ested in exploring the effects of turbulence on macroscopic bodies 
of sizes > 100 m, and so it is useful for our discussion to assume 
an optimistic outcome for this model. Toward the end of this sec- 
tion we also discuss some of the shortcomings of our model, and 
issues that are raised by these for the interpretation of our results. 



7.1 Turbulent discs without dead zones 

The results presented in Section |5] suggest that planetesimal for- 
mation and growth via a process involving mutual collisions be- 
tween smaller bodies is not possible in fully turbulent discs. The 
rapid excitation of large random velocities which exceed the dis- 
ruption/erosion threshold for planetesimals with ^ 100 km will 
simply lead to the destruction of bodies which form in this manner. 
Although the turbulence simulated in this paper is quite vigorous. 



because of the imposed vertical magnetic field, the scaling devel- 
oped in paper I for the strength of stochastic forcing as a function of 
the effective viscous stress suggest that even an order of magnitude 
decrease in the effective viscous a value will not decrease the ran- 
dom velocities sufficiently to prevent catastrophic disruption from 
occurring. 

The formation of large {R^ =^ 100 km) planetesimals through 
chondrules concentrating in turbulent eddies may be possible in 
fully turbulent discs driven by the MRI. The obvious requirement 
for there to be a turbulent cascade resulting in Kolmogorov- scale 
eddies in which chondrules can concentrate would appear to be 
satisfied in such a disc. During the initial formation and settling 
stage of these objects, they are likely to be of low density and hence 
strongly coupled to the gas via gas drag. Relative velocities on lo- 
cal scales relevant for collisions are likely to be small. As they con- 
tract to form "sandpile" planetesimals, however, they will decouple 
from the gas and evolve dynamically like planetesimals with inter- 
nal densities of pp ^ 2 g cm"' . If formation/contraction times for 
these bodies are ^ 10** yr at 5 au, then turbulent stirring will cause 
their velocity dispersion to grow to the surface escape velocity of 
=^ 100 ms"' within 3.5 x 10''yr, preventing runaway growth from 
ensuing. Further collisional growth between sandpile planetesimals 
will therefore occur slowly. Continued driving of the velocity dis- 
persion by turbulence eventually allows v^jsp to exceed the disrup- 
tion value for 100 km rubble piles within approximately 1.2Myr. 
Clearly collisional growth must allow significantly larger bodies to 
form within this time to prevent the eventual collisional destruction 
of the majority of sandpile planetesimals. Even if larger bodies do 
form that are safe against collisional disruption, questions remain 
about the possibility of forming large planetary cores that can ac- 
crete gas within a few Myr in the absence of runaway growth. A 
model in which a relatively small number of large planetesimals 
avoid collisional destruction and grow by accreting the surround- 
ing collisional debris may, however, provide a viable route to the 
growth of massive planetary cores. 

Rapid formation of large planetesimals through the clump- 
ing of metre-size boulders, followed by their gravitational collapse, 
has been demonstrated in f ully turbulent discs driven by the MRI 
djohansen et al.l2007l . l20T l|). Although the size distribution of plan- 
etesimals arising from this process is not known accurately because 
of numerical limitations, the simulations produce objects with sizes 
similar to Ceres (^p ^ 500 km). Such objects have surface escape 
velocities of 500 ms"', and the time over which turbulent stirring 
generates a velocity dispersion of this magnitude is ~ 8 x 10^^ yr. It 
is unclear whether a population of planetesimals with an approxi- 
mately unimodal size distribution centred on Rp = 500 km can un- 
dergo runaway growth, because of self-stirring and weak gas drag 
damping of eccentricities. But if it is feasible then turbulent stirring 
in a fully active disc such as computed in model Al will probably 
not provide significant hindrance because of the long time scales re- 
quired for the velocity dispersion to grow above v^isp =^ 500 m s"' . 

7.2 Tlirbulent discs witli dead zones 

Our results for model Dl show that in a model with a relatively 
deep dead zone, the equilibrium velocity dispersion is determined 
by collisional rather than gas drag damping for bodies with R^, < 
10 km, for reasons discussed in section [53] The equilibrium ve- 
locity dispersion lies between the disruption/erosion thresholds for 
weak and strong aggregates, suggesting that incremental collisional 
growth of planetesimals is possible, but depends on the material 
strength of the bodies. For 100 m bodies, the time required to ex- 



© 2002 RAS, MNRAS OOO.fTlfTsl 



Dynamics of planetesimals in dead zones 15 



cite the velocity dispersion to thie equilibrium value of 4.5 ms"' is 
2.2 X 10'* yr, comfortably shorter than the formation time of km- 
sized planetesimals at 5 au. Similarly, the time required to excite 
1 km bodies to their equilibrium random velocities is 2.2 x 10^ yr, 
comparable to our assumed formation time for these bodies through 
incremental growth. At each stage during the incremental growth 
of planetesimals within the dead zone, the planetesimals will have 
velocity dispersions close to the equilibrium values, which then 
control in-part the formation and growth rates. But provided the 
material strength of the planetesimals is sufficient (larger than that 
for weak aggregates), planetesimal growth should be possible in a 
dead zone through coUisional agglomeration, s o long as the bounc- 
ing and metre- si ze barriers can be overcome dSrauer et al.11200^ ; 
IZsometalJl2010l) . 

Runaway growth of km-sized planetesimals requires the ve- 
locity dispersion to be significantly lower than the escape veloc- 
ity from the accreting bodies, and for 10 km planetesimals Vesc = 
10ms"'. Results from model Dl indicate that a velocity disper- 
sion of this magnitude will be excited within ^10^ yr, effectively 
quenching the runaway growth. Continued collisional growth can 
still proceed, since the velocity dispersion remains below the dis- 
ruption/erosion threshold (at least for strong aggregates and rubble 
piles), and the equilibrium velocity dispersion for 10 km objects 
I'disp = 51ms"' is not reached until after 2.2 Myr. If growth ensues 
quickly enough, the eventual formation of bodies with Rp- 100 km 
will allow runway growth of these planetesimals to occur. Model 
Dl shows that the time required to raise the velocity dispersion 
up to the escape velocity I'esc = 100ms"' for these bodies is 
10 Myr. In this general picture, the influence of stochastic gravi- 
tational forces acting on planetesimals in a dead zone can cause the 
postponement of runaway growth for 10 km planetesimals, but once 
100 km sized bodies form, rapid runaway growth produce numer- 
ous ~ 1000 km sized oligarchs which will then undergo oligarchic 
growth to form planetary embryos and cores. This general picture 
of delayed runaway growth is similar in concept to that suggested 
bv lKortenkamp et al ( 1200 ih for planet formation in binary systems. 

The formation of sandpile planetesimals via the concentra- 
tion of chondmles within turbulent eddies should be possible in 
a dead zone, provided that the Reynolds stress generated there 
is indicative of a turbulent flow with an energy cascade resulting 
in Kolmogorov-eddies in which millimetre particles can become 
trapped. At the present time it is not clear that this is an accurate 
description of the dead zone flow, as the Reynolds stress appears 
to be largely due to spiral density waves that are stochastically ex- 
cited in the disc active regions, and undergo non linear damping 
as th eir wavelengths shorten during ra dial propagation through the 
disc l lHeinemann & PaDaloizoti2009bl ). A detailed analysis of high 
resolution simulations is required to examine the nature of the flow 
in the dead zone. Assuming that 100 km sandpile planetesimals can 
form, the slow growth of the velocity dispersion means that turbu- 
lent stirring cannot provide an impediment to subsequent runaway 
growth. 

Similar comments apply to the rapid formati on of planetesi- 
mals i n dead zones via the mechanism described inljo hansen et al.l 
( |2007|) . Formation in this case is driven to a large degree by trap- 
ping of metre-size bodies in pressure maxima (lo calised regions 
of an ticyclonic vorticity) generated by turbulence jjohansen et al.l 
l20Tlh . and these are not as prominent in the midplane of a disc 
with a dead zone. This suggests that the stream ing instability mod- 
els in laminar discs presented in ljohansen et al. k2009) may provide 
a more promising avenue for rapid planetesimal formation in dead 
zones. Once large planetesimals form, the possible onset of run- 



away growth will not be inhibited by the weak stochastic excitation 
of the velocity dispersion in a dead zone. 

7.3 Caveats 

We now consider the possible effects of assumptions we have made 
in our models that may affect their outcomes and the conclusions 
we have drawn from them. 

7.3.1 Disc mass 

The disc model we adopt in this study is lower in mass by a factor 
of ~ 3 than disc models that have been used as the basis for suc- 
cessful models of giant planet formation jPoUacket all 19961) . But 
it should be noted that these models often neglected the effects of 
planetary migration which can enable planetary embryos to grow 
beyond their isolation masses, thus assisting in the formation o f 
giant planets within expected disc life times jAlibert et al.ll200^ . 
If relative density perturbations in the dead zone remain the same 
when the disc mass is increased, then a higher disc mass will lead 
to a corresponding increase in turbulent stirring, making planetes- 
imal destruction/erosion more likely. But, in a disc where the ion- 
isation sources are external (stellar X-rays, galactic cosmic rays), 
the column density and mass of the active zone are almost indepen- 
dent of the total disc mass/surface density, and depend primarily 
on the penetration depth of the ionising sources. Consequently, we 
may expect the relative density fluctuations near the midplane in the 
dead zone to be weaker in a heavier disc than in a lighter disc, po- 
tentially providing a less destructive source of gravitational stirring. 
At present it is unclear how turbulent stirring in a dead zone scales 
with the disc mass, and studying this issue is difficult from the com- 
putational point of view. The shearing box models we present here 
are extended in the x and y directions compared to those used in 
most studies, and include ±5.5 scale heights above and below the 
midplane. Increasing the disc mass moves the boundary between 
the dead and active zones further from the midplane, such that the 
vertical computational domain must be extended significantly to ac- 
commodate the active zone, increasing the computational expense. 
It is our intention to address this issue in a future study. 

7.3.2 Initial magnetic fields 

We have assumed a particular configuration for the initial magnetic 
fields in our simulations that includes a small but non-negligible 
net-flux vertical component. It is well known that the existence 
and strength of net-flux magnetic fields changes the saturation 
level of developed MHD turbulence in simulations of the MRI 
( lHawlevetal]|l995l) . A future study of planetesimal dynamics em- 
bedded in turbulent discs with dead zones should examine ffie role 
played by the field topology and strength. 

7.3.3 Planetesimal size distribution 

In our estimates of collisional damping we assume that all plan- 
etesimals are the same size, and that all solids are incorporated into 
planetesimals. The latter assumption is broadly consistent with our 
disc model, in which we assume that 90% of the dust has been 
depleted due to the growth of grains into larger bodies. But the 
assumption of a single planetesimal size is not realistic. A more ac- 
curate calculation of the role of collisional damping would require 
a self-consistent model of accretion and fragmentation such that the 
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size distribution is modelled self-consistently. Unfortunately such 
a model goes beyond the scope of this paper. 

7.3.4 Planetesimal scale height 

When estimating the scale height of the planetesimal disc dur- 
ing the calculation of collisional damping rates, we assumed that 
Hp = Vdisp/f2k (i-e. the velocity dispersion is isotropic). In a turbu- 
lent disc, however, where a significant component of the stochastic 
gravitational field is generated by spiral density waves that have lit- 
tle vertical structure, we should expect that the eccentricity growth 
exceeds that of the inclination. Our simulation results, however, 
show that for models Dl and D2 the inclination growth actually 
exceeds that of the eccentricity. This suiprising result arises be- 
cause the disc midplane (centre of mass in the z-direction), oscil- 
lates vertically with an amplitude that reaches ~ 0.02/f , causing an 
oscillation of similar magnitude in the planetesimals. Future simu- 
lations of planetesimals within dead zones will be computed using 
global models to examine how robust this phenomenon is, although 
we note that a similar vertical oscillat ion has been repor ted in the 
shearing box simulations presented bv lOishi et al.l ( |2007|) . 

7.4 Equation of state 

Wave propagation in discs is strongly influenced b y the equation of 
state jLin et al.ll990l : lKorvcanskv & Pringlell995l) . This may affect 
the ability of waves that are excited in the active regions to propa- 
gate into the dead zone, and thus affect the turbulent stirring there. 
The simulations presented in this paper all adopted an isothermal 
equation of state, whereas protoplanetary discs are expected to be 
optically thick at 5 au. 

For parameters typical of T Tauri stars, the primary heating 
mechanism near the midplane is the absorption of the radiation 
emitted by the disc atmosphere, whi ch in turn heats mainly b y 
absorbing light from the central star dChiang & Goldreichlll997h . 
The accretion heating also occurs in the atmosphere, leading to 
a thoroughly isotherm al interior when a dead zone is present 
( iHirose & Tumeilboilh . The role of the equation of state in de- 
termining the level of turbulent stirring at the midplane will be ex- 
plored in a future publication. 



8 CONCLUSIONS 

In this paper we have presented the results of simulations of plan- 
etesimals embedded in local models of protoplanetary discs. The 
main aim of this work is to understand the dynamics of planetesi- 
mals in turbulent discs, and to examine how dead zones of varying 
vertical thicknesses change the stochastic forcing of the radial ve- 
locity dispersion, and the diffusion of planetesimal semimajor axes. 
We consider the influence of gas drag and collisional damping in 
counteracting the growth of random velocities due to turbulent stir- 
ring, and estimate the equilibrium velocity dispersions which re- 
sult. The main conclusions of this study are as follows. 

The presence and depth of a dead zone have negligible effect 
on the measured coherence time of the particle stirring, but only 
affect its amplitude as measured through the width of the torque 
distribution. 

In fully turbulent discs without dead zones, we find that the 
velocity dispersion of embedded planetesimals grows rapidly, and 
quickly exceeds the threshold for disruption for planetesimals in the 



size range 100 m < Rp < 10 km. We conclude that planetesimal for- 
mation via collisional accretion of smaller bodies cannot occur in 
globally turbulent discs. The direct formation of large 100 km-sized 
planet esimals via turbulent concentration and gravitational insta- 
bility jjohansen et al.ll2007l ; ICuzzi et al]|2008l) may occur, but the 
excitation of their random velocities may prevent runaway growth, 
making it difficult to form large planetary embryos and cores within 
reasonable disc life times. 

Radial diffusion of planetesimals occurs over distances of 
~ 2.5 au in fully turbulent discs over typical disc life times of 
5 Myr. This is inconsistent with the observed correlation between 
asteroid taxonomic types and heliocentric distance in the asteroid 
belt, which indicate radial mixing over distances < 0.5 au during 
the formation of the solar system. It thus seems implausible that 
the solar nebula sustained levels of turbulence observed in numeri- 
cal simulations of fully turbulent discs. 

The stochastic forcing of planetesimal random motions is 
much weaker in discs with dead zones. In a dead zone with ver- 
tical extent ±2H above and below the midplane, the stochastic ex- 
citation of planetesimal eccentricities is reduced by a factor ~ 17 
relative to a fully active disc. The damping of planetesimal ran- 
dom motions in this model is dominated by collisional damping, 
and the equilibrium velocity dispersions lie between the disruption 
thresholds for weak and strong aggregates. It thus appears that the 
collisional growth and formation of km-sized planetesimals may 
be possible in a dead zone . But, we note that the tidal break-up of 
comet Shoemaker-Levy 9 jAsphaug & Benzll996l) . and analysis of 
the colli sion ejecta from comet Tempel 1 during the Deep Impact 
mission jRichardson et alj2007l) . suggest that cometary bodies pos- 
sess modest material strength. It therefore remains uncertain what 
range of collision velocities can be tolerated during the incremen- 
tal growth of km - sized plan etesimals as envisag ed in the models of 
IWeidenschiliin3 ( |2000|) and lBrauer et"ai] ( |2008[) . 

Although km-sized planetesimals may be able form within a 
dead zone, subject to overcoming the metre-sized and bouncing 
baniers, stochastic forcing by the turbulent disc excites random ve- 
locities above the escape velocity for = 10 km bodies within 
10'' yr, quenching runaway growth. Slower growth of planetesimals 
through collisions remains possible, and runaway growth can ensue 
once bodies with Rp - 100 km form, due to the slow growth of ran- 
dom velocities within a dead zone. If Rp ^ 1000 km 'oligarchs' 
form, oligarchic growth can proceed with very modest influence 
from the weak stochastic forcing in the dead zone. A corollary of 
this is that 100 km-sized planetesimals that form through the con- 
centration of particles and subsequent gravitational instability in a 
dead zone will not be prevented from undergoing runaway growth 
by turbulent stirring. 

In summary, we come to the following specific conclusions 
regarding the outcome of our present study: 

• Comparing the results obtained for dead zones with different 
vertical heights (~ 1.5// and 2H), we find that the larger dead zone 
model clearly favours the collisional formation of planetesimals be- 
cause of weaker stochastic forcing of random motions. 

• The radial diffusion of planetesimals is much reduced in discs 
with dead zones relative to fully active discs. Diffusion over a dis- 
tance equal to 0.26 au is predicted by our deeper dead zone model, 
consistent with solar system constraints. 

• The low diffusion levels for planetesimal semimajor axes in 
dead zones clearly means that stochastic torques in such an envi- 
ronment cannot prevent the large scale inward drift of low mass 
planets via type I migration. 
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• We suggest that dead zones provide safe havens for km-sized 
planetesimals against destructive collisions, and are probably a re- 
quired ingredient for the formation of planetary systems. 

This paper is part of series in which we are examining the 
dynamics of planetesimals and planets in increasingly realistic disc 
models, and shows that different qualitative outcomes for planet 
formation can be expected in discs with dead zones versus those 
without. In future publications we will examine a number of issues, 
including how changing the disc surface density and the mass ratio 
between the active and dead zones modifies the turbulent stirring, 
and how turbulence varies as a function of heliocentric distance in 
global MHD simulations of protoplanetary discs. 
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